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1. IWlMXJUCriCN 


Although from a historical point of view, composite materials have 
found practical use for centuries, during the past two decades there has 
been a trenendous increase in the use of composite materials in engineering 
applications, particularly in aerospace engineering. One of the biggest 
attractions of these materials is that one is able to design and 
manufacture such materials to sustain a specific type of loading in a most 
efficient manner. If properly produced, composites can often achieve a 
combination of properties that are far superior to the properties of the 
individual constituents acting independently. 

It is well known that gas turbine engine structures, particularly 
those components directly in the hot gas flow path, are subjected to 
extremely severe thermal and mechanical loading that can often lead to 
creep enhanced distortion, cracking and low cycle fatigue. As the demand 
for more efficient propulsion system rises so does the thermal and 
mechanical loading. It is unlikely that the current generations of metal 
alloys would be suitable candidates for structural components in the future 
generations of efficient propulsion systems. Ceramic components are often 
thought to be ideal as far as their thermal durability is concerned. 
Unfortunately, ceramics do not have adequate tensile strength to sustain a 
high level of mechanical loading. In recent years there has been 
significant effort in the attempt to incorporate fibrous inclusions within 
a ceramic natrix to develop a class of new materials (ceramic composites) 

for advanced engineering applications. 

The mechanical behavior of ceramic composites under nonlinear, thermal 
and dynamic loading is extremely complex and can only be understood if the 
observed behavior is interpreted in terms of micromechanical analyses. 
Such analyses must take care of the complex interaction of the individual 
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fibers or bundles of fibers embedded in the three-dimensional ceramic 
matrix and must allow for increasing levels of sophistication in terms of 
the idealization of the fibers as well as the ceramic matrix. In addition 
complex interface behavior and controlled failure of the fiber must be 
considered. 

This report details progress made during the period of March 1988 to 
December 1988 in a five-year program commencing in March 1988, towards the 
development of a boundary element code designed for the micromechanical 
studies of advanced ceramic composites. Additional effort has been made in 
generalizing the implementation to allow the program to be applicable to 
real problems in the aerospace industry. 

The primary goal of the first year has been to develop the boundary 
integral formulation for a fully-bonded elastic inclusion within an elastic 
matrix, and to implement the formulation in a boundary element program in a 
sufficiently general manner as to facilitate implementation of future 
development for the remainder of the present program of research. 

Significant progress has been achieved during the first year of the 
present effort. The analytic and numerical basis for the ceramic composite 
program has been developed. Ihis effort included: 

1. the derivation of the boundary integral formulation modified for holes 
and inserts, 

2. the derivation of the kernel function for one-dimensional line 
integration of holes and inserts, 

3. the implementation for the assembly and solution of inserts, and 

4. the validation and verification runs using the developed computer 
code. 
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The ceramic composite formulation has been implemented in the three- 
dimensional boundary element computer code 'BEST3D' which was developed for 
NASA by Pratt and Whitney and the State University of New York at Buffalo 
under contract NAS3-236 97. BEST3D has been adopted as the base for the 
ceramic composite program, so that many of the enhanced features of this 
general purpose boundary element code can be utilized. Some of these 
facilities include sophisticated numerical integration, the capability of 
local definition of boundary conditions, and the use of quadratic shape 
functions for modeling geometry and field variables on the boundary. The 
multi-region implementation permits a body to be modeled in substructural 
parts; thus dramatically reducing the cost of the analysis. Furthermore, 
it allows a body consisting of regions of different ceramic matrices and 
inserts to be studied. 

In the next section, the existing approaches for the study of the 
micromechanical behavior of composites are briefly reviewed. This is 
followed by the development of the boundary element formulation for ceramic 
composites in Section 3. Section 4 outlines the development of the general 
computer program implementation followed by the description of the 
program's data input and output in Section 5. Example data and results are 
included to demonstrate the convenience in modeling and analyzing 
composites using this code. In Section 6, a number of numerical examples 
are presented to demonstrate the power of the present implementation. This 
report is then concluded with a summary and plan for future development. 
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2. REVIEW OF EXISTING APPROACHES 

A number of techniques of varying degree of sophistication are 
presently used to study the micromechanical behavior of composites. In 
essence the various micromechanical analyses recognize the inhomogeneous 
nature of the composite material, but generally ignore finer details of the 
structure of the fiber and matrix to a varying degree. Numerous 
approximations are made concerning the packing geometry and response fields 
within the body so that mathematical analyses can be performed to relate 
the properties and concentration of fibers and matrix to the average 
property of the body. A complete bibliography of these methods which are 
variously called 'mechanics of materials', 'self-consistent fields', 
'variational methods' and ’numerical techniques' are given in Chamis and 
Sendeckyj (1968). 

Although none of the above mentioned work specifically focuses on 
ceramic composites (and it appears that very little is available in the 
current literature), the overall conclusions derived apply reasonably well 
to the present case. 

Rational methods based on mechanics of materials have been developed 
over the last decade to explain and predict the behavior of composite 
systems in terms of their materials make-up. Some of these simple models 
have indeed been very effective. For example, a simple model of the fiber 
and resin components responding by parallel reactions to imposed thermal or 
mechanical loads is generally adequate to describe and predict the 
following longitudinal properties of unidirectional composites: 

- Longitudinal modulus E* J 

- Poisson's ratio, 

- Linear thermal expansion coefficient, a^; 

- Thermal conductivity, K L * 



Primary emphasis in the model has been given to stiffness limited 
structures so that attention has been focused on the longitudinal modulus. 
For all practical purposes the contributions of the matrix phase to this 
property are often assumed negligible. 

The properties perpendicular to the direction of the fibers are not so 
simply described nor so readily predicted. A simple model of the fiber and 
matrix components responding in series to imposed loads leads to 
conservative estimates of the properties in this direction. More 
comprehensive and elaborate models show that the transverse properties and 
shear moduli are sensitive to: 

- the shape of the fibers; 

- the packing geometry of the fibers; 

- the variations in spacings of the fibers; 

- the properties of the matrix. 

Ambiguities in the theoretical models make it difficult to isolate and 
identify the extent to which each of these factors influences the 
transverse properties and shear moduli of composites; qualitatively, it is 
clear that the matrix dominates the behavior of these properties. 

A model for compressive strength based on a failure mode of in-phase 
fiber buckling leads to the conclusion that the compressive strength of 
unidirectional composites is determined by the shear modulus of the matrix. 
Experimental observations can be rationalized in terms of this model if the 
effective shear modulus of the matrix is assumed to be less than that of 
the comparable bulk material. Adherence to this model further requires 
that the effective shear modulus of the matrix be dependent upon the system 
of reinforcing fibers used and/or dependent upon the degree of adhesion 
between the fibers and the matrix. 
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Current theoretical models cannot unambiguously explain or predict the 
tensile strength of unidirectional composites; however, examination of 
models for the extremes of fiber-matrix coupling leads to a qualitative 
identification of the critical factors that influence tensile strength: 

- the statistics of the fiber strengths; 

- load transfer efficiency; 

- resistance of the natrix to crack propagation; 

- the adhesive bond strengths between the fiber and the matrix 

The direct contribution of the strength of the natrix to the tensile 
strength of the composite is negligible, nonetheless, an imperfect matrix 
can seriously detract from the realization of the full strength potential 
of the fibers through indirect influences involved in load-transfer 
efficiency or crack sensitivity. 

Load-transfer efficiency and crack sensitivity appear to be diametric 
functions of the adhesive bond strength between the fiber and the matrix. 
This implies that some optimum (not necessarily a maximum) adhesive bond 
strength is required to establish a proper balance between load transfer 
efficiency and the overall strength-toughness of the composite. 

It was suggested in some physiochemical models that the molecular 
mechanisms involved in achieving adhesion between the fiber and matrix 
could significantly perturb the molecular structure of the matrix; these 
structural perturbations may develop in an interphase region whose 
properties differ appreciably from the properties of bulk material. 
Additionally, stress diffusion from fiber to the matrix may also alter the 
properties of the matrix locally. Since very large surface areas are 
involved in the contact between the fiber and the matrix, a sufficient 
quantity of interphase material could be generated to influence the average 
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in situ properties of the matrix. Hence, it is reasonable to expect that 
the average properties of matrix in composites may differ from the 
corresponding properties of bulk material and may differ between various 
reinforcing fiber systems. 

In conclusion, the contributions of the reinforcing fibers to the 
performance of composite materials are qualitatively well understood. 
Fiber dominated properties (i.e., the longitudinal modulus) can be 
adequately explained and predicted. The direct contribution of the matrix 
to the longitudinal modulus and tensile strength is negligible; however, 
secondary effects associated with the matrix can detract frcm the potential 
tensile strength of the system. Although current simple theoretical 
models indicate that the matrix controls the compressive strength, 
transverse moduli, and shear moduli of composites, this is not sufficient 
for the solution of boundary value problems involving nonlinear, thermal, 
and dynamic loading. 

Comparatively little effort has been given to identify the specific 
roles of the matrix in the overall performance of advanced composite 
materials. This is largely due to the emphasis cm longitudinal modulus as 
well as the theoretical and experimental difficulties involved in isolating 
the contributions of the matrix from effects due to fiber shape and packing 
geometry. The possibility that the properties of the matrix can differ 
from the properties of bulk material, as well as differ between various 
fiber systems, further complicates the in situ matrix. 

Since the interaction of the fiber and matrix is complex, it is likely 
that the best technique to use to study the micromechanical behavior would 
be a numerical method. It is therefore not surprising that the finite 
element method is often used to develop the micromechanical model. 
Although it is perhaps the most powerful numerical method over the entire 
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spectrum of engineering science, it has never been fully satisfactory for 
problems with high stress gradient, complex interface phenomena, high 
thermal gradients and large variations in the stiffness within the same 
body. Unfortunately all of these occur within a ceramic composite 
assembly. Nevertheless, the use of finite element to study the 
micromechanics of composites is quite common, perhaps primarily due to the 
absence of a reasonable alternative. 
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3. BOCHiAK? ELEMENT FORMULATION FOR CERAMIC COMPOSITES 

3 . 1 Introduction 

It is evident that for proper micromechanical analysis of ceramic 
composites one needs to use a numerical method that is capable of 
idealizing the individual fibers or individual bundles of fibers embedded 
within a three-dimensional ceramic matrix. The analysis must be able to 
take account of high stress gradients from diffusion of stress from the 
fiber to the ceramic matrix and allow for the interaction between the 
fibers through the ceramic matrix. The analysis must be sophisticated 
enough to deal with failure of fibers described by a series of increasingly 
sophisticated constitutive models. Finally, the analysis must deal with 
micromechanical modeling of the composite under nonlinear thermal and 
dynamic loading. 

The boundary element method is uniquely suited for the task. BEM has 
proven its ability to accurately determine stress near stress 
concentration. All functional quantities in a BEM system are on the 
boundary and interface surfaces, therefore, allowing nonlinear interaction 
between the matrix and insert interface to be readily described by failure 
models. Furthermore, recent development has shown the generality and 
versatility of boundary element method in analyzing large two- and three- 
dimensional models subjected to static, dynamic and thermal loads involving 
materials with nonlinear behavior. 

3.2 Boundary Integral Equations 

The conventional boundary integral equation for displacement is the 
starting point for the ceramic composite formulation. Ihe displacement for 
a point £ inside the elastic composite matrix is given below 

C.j(5) Ui (?) = j [G i j(x,?)t i (x) - FjjCx.SJu^xndStx) 

S 


9 



N 



where 




are the fundamental solutions of the governing differential 
equations of the ceramic natrix of infinite extent 


C are constants determined by the relative smoothness at 5 
u^.t^ are displacements and tractions 

S, s n are surfaces of the matrix and holes (left for fiber) 
respectively 


N is the number of fibers 


The conventional boundary integral equation for displacement can also 
be written for each of the N insert fibers. For the displacement at a 
point l inside the m^ insert we can write 

C™U)G.U) = J [Cjf j (x.5)t i (x) - E? j (x.l)u i (x)]dS m (x) (3.2) 

i.j = 1.2,3 

nP f 1 ? • are the fundamental solutions of the m^ 1 insert 

ij i] 

q. are constants determined by the relative snoothness at 5 in 

insert m 

u-.L are displacement and tractions associated with the m^ 
insert 

S* the surface of the m^ 1 insert 
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We next examine the interface conditions between the composite matrix and 
the insert. For a perfect bond the displacement of the matrix and the 
displacement of the inserts along the interface the interface are equal and 
the tractions are equal and opposite. 


Uj(x) = Uj (x) (3.3a) 

tj(x) = - tj(x) (3.3b) 

For a stiff insert in which the elastic modulus is much greater than the 
modulus of the composite matrix, the Poisson ratio of the insert can be 
assumed equal to that of the matrix with little error. Therefore, upon 
consideration of the surface normals at the interface and examination of 
the F^ kernels, we can write the following relation for the m insert 

^(x,?) = - F^j(x,$) (3 - 3c) 

Substitution of equations (3.3) into equation (3.2) yields the following 
modified boundary integral equation for insert m. 

C^.(?) Ui («) - J s m l - G ij (x ’5 )t i ( ^ ) + F i j(x.5)u i (x)]dS m (x) (3.4) 


Finally adding N insert equations (3.4) to equation (3.1) and cancelling 
terms, yields the modified boundary integral equation for the composite 
matrix 




= J lG i j(x,5)t i (x) - F i j(x,Uu i (x)ldS(x) 


N 

+5 J tG i j(x,?)t i (x)dS n (x) 
n=l * 

where 


(3.5) 


- Gij (x,?> - <$,<*.«> 


11 



C.j constants dependent on the geometry at ? 

3.3 Nunerical Implementation 

3.3.1 Discretization 

The integral representations of the previous section are exact 
statements of the ceramic composite problem, however, approximations such 
as finite discretization and numerical integration are necessary in order 
to obtain a solution to non-trivial problems. The goal of the numerical 
implementation of the present formulation is to obtain the most accurate 
and efficient implementation possible. 

The first step in this process is the conversion of the two- 
dimensional kernel integral over the surface of the hole and insert into a 
line integral. By performing an analytical integration in the 
circumferential direction on the surface of the hole (or insert) a 
considerable amount of computational time can be saved in the numerical 
integration. In this process the holes and inserts are assumed to be 
circular and a circumferential variation of a o +ajcos0 + a 2 sine is assumed 
in the displacements and tractions on the surface of the holes and inserts. 
Furthermore, tensor transformations on the kernels are necessary for 
inserts oriented at oblique angles with respect to the global axes. The 
resulting kernels, which contain a large family of elliptical integrals, 
are long and formidable. The complexity of these kernels prohibits their 
presentation in a tidy manner, and therefore, they have not been presented 
in this report. 

Once the analytical integration is complete, the inserts are 
discretized in their axial direction using linear or quadratic shape 
functions. In discretized form, equation (3.4) can be written for a single 
insert as 
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(3.6) 


P 

<H-j<«> u i (? > " ' 2 [ L^ j <*.5» )T( ’i> asP 3 tj 

p-i 

+ l [ I p F ij <*.5»J T <n>asP 3 “I 

P-1 s 

where P is the number of line elements, and 

N^n) represents the shape function across the line element. 
Summation over r is implied. 

In a similar rranner, equation (3.5) can be discretized using one- and 
two-dimensional shape functions in the following manner. 

Q 

Cij( ? )u i ( 5) = } [ } q G. j (x.?)M P (ti 1 . 1l2 )dS q ] t[ 

q=i s 

" 1 [ j q F ij (x.5)M p (ii 1 .n 2 )dS q ] uf (3.7) 

q=l S 

+ 1 [ L p G ij (x,ON Y (Ti)dsP ] tT 
p=l 

where Q is the number of surface elements on the outer surface of the 
composite matrix of the region, and 

(t| t^) represents the two- dimensional shape function. 
Summation over y and p is implied. 

Note that the same number of nodes, and consequently shape functions, 
are not necessarily used to describe both the geometric and functional 
variations. Specifically, in the present work, the geometry is exclusively 
defined by quadratic shape functions. On the other hand, the variation of 
the primary quantities can be described, within an element, by either 
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quadratic or linear shape functions. (The introduction of linear 
variations provides computationally advantageous in some instances.) 

3.3.2 tftmerical Integration 

The complexity of the integral in the discretize equation necessitates 
the use of numerical integration for their evaluation. The steps in the 
integration process for a given element is outlined belcw: 

1. Using appropriate Jacobian transformations, the curvilinear line 
elements and surface elements are mapped onto a unit line and planar unit 
cells, respectively. 

2. Depending on the proximity between the field point (?) and the 
element under consideration, there may be element subdivision and 

additional mapping for improved accuracy. 

3. Gaussian quadrature formulas are employed for the evaluation of 
the discretized integral over each element (or sub-element). These 
formulas approximate the integral as a sum of weighted function values at 
designated points. The error in the approximation is dependent on the 
order of the (Gauss) points employed in the formula. To minimize error 
while at the same time maintaining computational efficiency, optimization 
schemes are used to choose the best number of points for a particular field 
point and element (Watson, 1979). 

4. When the field point coincides with a node of the element being 
integrated, the integration becomes singular. In this case, the value of 
the coefficients of the F.. kernel corresponding to the singular node 
cannot be calculated accurately by numerical integration. Instead, after 
the integration of all elements is complete, this value is determined so as 
to satisfy a rigid body displacement of the body (Banerjee and Butterfield, 
1981). 
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3.3.3 Assembly of Equations for Composite Inserts 

After the derivation of the modified boundary integral equations and 
the analytical circumferential integration of the kernel functions, the 
next critical step in the formulation is the assembly of the inserts into 
the system equations. Here, efficiency is of utmost importance. The 
approach to writing an efficient algorithm is to keep the number of system 
equations to a minimum by eliminating all unnecessary unknowns from the 
system. The strategy used is to retain in the system only traction 
variables on the matrix/ insert interface. This is in contrast to a general 
multi-region problem where both displacement and tractions are retained on 
an interface. The elimination of the displacements on the interface is 
achieved through a backsubstitution of the insert equations into the system 
equations which are made up exclusively from equations written for the 
composite matrix (on the outer surface and on the surface of the holes). 
The procedure is described below. 

Equation (3.7) is used to generate a system of equations for nodes on 
the outer surface of the matrix and for nodes on the surface of the holes 
containing the inserts. Written in matrix we have 

On Matrix Surface: Gt - Fu + Gt? = 0 (3.8a) 

On Hole for Insert: Gt - Fu + Gt? = Iu? (3.8b) 

where 

t and u are traction and displacement vectors on the outer surface of 

the composite matrix 

t H and i? are traction and displacement vectors on the hole 

I is the identity matrix 

Our goal is to eliminate u^ from the system. To this end, equation (3.6) 
is written for each node on an insert, collocating slightly outside the 
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( 3 . 9 ) 


insert [where =0], we obtain 

Fu 3 = Gt 3 


Post multiplying equation (3.8b) by the F matrix in equation (3.9) yields 


FGt = FFu + FGt 3 « Fu 3 


( 3 . 10 ) 


Equation (3.9) can now be set equal to equation (3.10) and the final form of the 
system is derived. 


On Matrix: 
On Hole: 


Gt - Fu + Gt 3 =0 

FGt - FFu + (FG - G)t 3 = 0 (3.11) 


At every point on the outer surface, either the traction or the 
displacement is specified and on the surface of the hole only the tractions 
are retained. Therefore, the number of equations in the system are equal 
to the final number of unknowns, and hence, the system may be solved. 
Thereafter, equation (3.8b) is used to determine the displacement on the 
matrix/ insert interface. 

It should be noted that since the displacement about a particular hole 
is present only in the insert equation corresponding to that hole, 
backsubstitution can be performed one insert at a time in a more efficient 
manner than backsubstitution of all inserts at once. Further note that 
nowhere in the assembly process is a matrix inversion necessary. This 
efficient assembly process was made possible due to the unique formulation 
of the modified boundary integral equations developed earlier in this 
section. 

When the composite matrix is divided into a multi-region model, the 
above insert assembly is performed for each region independently. 
Thereafter, equilibrium and compatibility conditions are invoked at common 
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interfaces of the substructured matrix composite. After collecting 
together the known and unknown boundary quantities, the final system can be 
expressed as 

. Bhy (3.J2> 

where 

x is the vector of unknown variables at boundary and interface 
nodes, 

y is the vector of known variables, and 

A^.B^ are the coefficient matrices 

Standard numerical procedures can be used to solve the unknowns in 
equation (3.12). Details are described in the computer development 
section. 

3.4 Interior Quantities 

Once all the displacements and tractions are known on the matrix 
outersurface and on the matrix/ insert interface, interior quantities of 
displacement, stress and strain can be determined at any point in the 
composite matrix or in the insert. For displacement either the 
conventional boundary displacement integral equation (3.1) or (3.2) can be 
employed or alternatively the modified equations (3.3) or (3.5) can be 
used. 

Equations for strains can be derived from the forementioned 
displacement equations and the strain-displacement relations. Thereafter, 
equations for stress are obtained by substituting the resulting strain 
equations into Hooke's law. 

The resulting equations, however, are not only invalid on the surface, 
but also difficult to evaluate numerically at points close to it. For 
points on the surface, the stresses can be calculated by constructing a 
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local Cartesian coordinate system with the axes 1 and 2 directed along the 
tangential directions and the axis 3 in the direction of the outward 
normal. The stresses referred to these local axes (indicated by 

overbars) are then given by: 

°li = if^ ^3 + ^2 (e ll + e 22 5 + l+"v e u 
- _ E - 

a 12 “ a 21 - 2(l+v) s 12 


°22 = if ; l 3 + *_2 (8 11 + C 22 } 


'11 


E - 
1+ v e 22 



(3.13) 


a 32 " °23 ” ^2 
a 31 = a 13 = *1 

where E is the Young's modulus, e^j defines the components of the strains 
in the local axes system and t i are the traction on the boundary. This 
method of evaluating the stresses on the surface was originally devised by 
(Rizzo and Shippy, 1968). 
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4. COMPUTER PROGRAM DEVELOPMENT 

4.1 Introduction 

The goal of the computer program developed for ceramic composites is 
the accurate and efficient implementation of the formulation described in 
Section 3. Of equal importance is the degree of generality required in the 
definition of component geometry, loading and material properties. This is 
necessary if the program is to be applicable to real problems in the 
aerospace industry. 

For this reason the ceramic composite formulation has been implemented 
in the three-dimensional boundary element computer code 'BEST3D' (Boundary 
Element Stress Technology - Three-dimensional) which was developed for NASA 
by Pratt and Whitney and SONY /Buffalo under contract NAS3-23697. Since its 
development, BEST3D has proven itself to be a highly accurate and 
numerically efficient boundary element program. 

The development of the computer program 'Ccmposite-BEST' is discussed 
in the following sections. 

4.2 Program Structure 

Composite-BEST is designed to be a fully general ceramic composite 
analysis system employing the boundary element method. The program is 
written using standard FORTRAN 77. Development has been carried out at 
SUNY/Buffalo on an HP9000 minicomputer system. The required code and 
workspace fit in core without requirement for overlays. The nature of the 
method is such that, for any realistic problem, not all required data can 
reside simultaneously in core. For this reason extensive use is made of 
both sequential and direct access scratch files. 

The program first executes an input segment. After the input has been 
processed, the surface integrals are calculated and assembled into the set 
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of system equations using specified boundary conditions, followed by the 
insert assembly and the inclusion of the insert equations into the general 
system. The system matrix is then decomposed and saved on disk, followed 
by the calculation of the solution vector. The full displacement and 
traction solution on each boundary element and insert element is then 
reconstructed from the solution vector. In a time dependent problem the 
process of constructing the load vector for the system equations is 
repeated at each time step, but the integration, formation and 
decomposition of the system matrix are done only once. 

Various aspects of the computer program are discussed below. 

4.3 Program Input 

The input for Composite - BEST is free field. Meaningful keywords are 
used to identify data types and to name particular data sets. The input is 
divided into five types: 

1. Case Control Cards 

The case control cards define global characteristics of the problem. 
In addition to the problem title, the times for multiple time steps are 
defined. The reading or writing of restart data is also defined at this 
point. The restart facility allows one to change the arrangements of 
fibers without recalculating the various coefficients. 

2. Material Property Definition 

The material property input allows the definition of material 
properties for a variety of materials. The Young's modulus can be 
prescribed in tabular form for a user-defined set of temperatures. 
Temperature independent values of Poisson's ratio are also defined. 
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3 . Geometry Input 

Geometry input is defined one GMR (generic modeling region, or 
subregion) at a time. To initiate the input, a tag is provided to identify 
the GMR, a material name and reference temperature are defined to allow 
initialization of material properties. 

The next block of geometry input consists of the Cartesian coordinates 
of the user input points for the outer surface geometry definition of the 
composite matrix, together with identifiers (normally positive integers) 
for these geometric nodes. 

Following the definition of an initial set of nodal points, the 
surface connectivity of the outer surface of the composite matrix is 
defined through the input of one or more named surfaces. Each surface is 
made up of a number of elements, with each element defined in terms of 
several geometric nodes. Three sided elements, defined using six rather 
than eight geometric nodes, are used for mesh transition purposes. The 
terms quadrilateral and triangle are normally used to refer to the eight 
and six noded elements, although the real geometry represented is, in 
general, a nonplanar surface patch. Seven and nine noded elements are made 
available by adding a central node to the six or eight noded elements. 

Over each element the variation of displacement and traction can be 
defined using either the linear or quadratic shape functions. Linear and 
quadratic elements can share a common side, which is then constrained to 
have linear displacement and traction variation. 

Finally an option is available to allcw quadratic functional variation 
(8 or 6 nodes) to be used in conjunction with linear geometry (4 or 3 
nodes). In this case the program generates the additional nodes 
automatically at mid-point of the sides. The characteristics of the 
various element types are summarized below. 
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Surface Element Type 


Geometry 

Nodes 


Displacement/Traction Nodes 


Linear Quadrilateral 8 or 9 

Linear Triangle 6 or 7 

Quadratic Quadrilateral 8 or 9 

Quadratic Triangle 6 or 7 

Quadratic Quadrilateral 4 

Gkiadratic Triangle 3 


4 

3 

8 or 9 
6 or 7 
8 

6 


Following the definition for the composite matrix outer surface, the 
embedded inserts are then defined. These are defined as curvilinear line 
elements with a prescribed radius of cross-section. The inserts are 
generally straight, however as noted, curved inserts are also allowed. The 
user first defines the nodal coordinates of the centerline of the insert. 
Thereafter, the radius and the insert connectivity is defined. Linear and 
quadratic elements are available for both geometry and functional 
variation, however, quadratic functional variation over linear geometry is 
not presently available. The various options for the insert elements are 
summarized be lew. 

Geometry 

Insert Element Type Nodes Displacement/Traction Nodes 

Linear-Linear 2 2 
Quadratic-Linear 3 2 
Quadratic-Quadratic 3 3 


Note only the surface of the insert needs to be defined, i.e., the hole in 
the composite matrix which encompasses the insert does not have to be 
explicitly defined. 


4. Interface Conditions 

The interface input describes the connection of surfaces or elements 
of one composite matrix region to another. Interfaces between the 
composite matrix and inserts ds not have to be defined. Special types of 
interface conditions which are available presently include fully-bonded and 
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sliding contact between two GMRs, and springs to other regions. In the 
current implementation, fully-bonded connections between the insert and 
matrix has been assumed. This will be relaxed in later work. 

5. Boundary Condition Input 

The final input section provides for the definition of boundary 
conditions, as functions of both position and time. Data can be input for 
an entire surface, or for a subset (elements or nodes) of a surface. Input 
can be in global coordinates, or can define rollers or pressure in the 
local coordinate system. Input simplifications are available for the 
frequently occurring cases of boundary data which is constant with respect 
to space and/or time variation. Each boundary condition set can be defined 
at a different set of times. 

4.4 Surface Integral Calculation 

Following the processing of the input data, the surface integrals 
occurring in equations 3.6 and 3.7 are evaluated numerically. This is the 
most time consuming portion of the analyses. In Compos ite-BEST the results 
of these integrations are stored as they are calculated, rather than being 
assembled into the final equation system immediately. Although this is 
somewhat more costly in terms of storage and CPU (central processing unit) 
time, it has led to much greater clarity in the writing of Composite-BEST. 
In addition, it provides much greater flexibility in the implementation of 
various restart and boundary condition options. 

The calculations proceed first by GMR (generic modeling region), then 
by source point (the equation being constructed) and finally by surface 
element and insert element. The results for each source point element pair 
are written to disk. All of the calculations are carried out and stored in 
the global (Cartesian) coordinate system. 
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Hie integration of the BEM equations is the most complex part of the 
code. In this process either singular or nonsingular integrals can be 
encountered. The integrals are singular if the source point for the 
equations being constructed lies on the element being integrated. 
Otherwise, the integrals are nonsingular, although numerical evaluation is 
still difficult if the source point and the element being integrated are 
close together. 

In both the singular and nonsingular cases Gaussian integration is 
used. The basic technique is developed in Banerjee and Butterfield, 1981. 
In the nonsingular case an approximate error estimate for the integral was 
developed based on the work of Stroud and Secrest (1966). Hiis allows the 
determination of element subdivisions and orders of Gaussian integration 
which will retain a consistent level of error throughout the structure. 
Numerical tests have shown that the use of 3, 4, and 5 point Gauss rules 
provide the best combination of accuracy and efficiency. In the present 
code the 4 point rule is used for nonsingular integration, and error is 
controlled through element subdivision. The origin of the element 
subdivision is taken to be the closest point to the source point on the 
element being integrated. 

If the source point is very close to the element being integrated, the 
nsp> of a uniform subdivision of the element can lead to excessive computing 
time. This frequently happens in the case of aerospace structures, due 
either to mesh transitions or to the analysis of thin walled structures. 
In order to improve efficiency, while retaining accuracy, a graded element 
subdivision was employed. Based on one-dimensional tests, it was found 
that the subelement divisions could be allowed to grow geometrically away 
from the origin of the element subdivision. Numerical tests on a complex 
three-dimensional problem have shewn that a mesh expansion factor as high 
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as 4.0 can be employed without significant degradation of accuracy. 

In each case of singular integration (source point on the elements 
being integrated) the element is first divided into subelements. The 
integration over each subelement is carried out using a Jacobian 
transformation in mapping. This coordinate transformation produces 
nonsingular behavior in all except one of the required integrals. Normal 
Gauss rules can then be employed. The remaining integral (that of the 
traction kernel times the isoparametric shape function which is 1.0 at 
the source point) is still singular, and cannot be numerically evaluated 
with reasonable efficiency and accuracy. Its calculation is carried out 
indirectly, using the fact that the stresses due to a rigid body 
translation are zero (Lachat and Watson, 1976). It has been found that 
subdivision in the circumferential direction of a two-dimensional surface 
element is required to preserve accuracy in the singular integration of the 
outer surface. A maximum included angle of 15 degrees is used. 
Subdivision in the radial direction has not been required. 

The integrals required for calculation of displacement and stress at 
interior or surface points are of the same type as those involved in the 
generation of the system equations, except that only nonsingular integrals 
are involved. If the source point involved is located on the surface of 
the body, then numerical integration is not required. Instead, the 
required quantities are calculated using the displacements and tractions on 
the element (or elements) containing the source point, as discussed in 
Section 3.4. 

4.5 System Matrix Assembly 

The first step in the assembly process is the reduction of the 
rectangular matrix of F integrals to a square matrix. This matrix is the 
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prototype of the system matrix. The columns of this matrix are transformed 
or replaced, as required by the boundary conditions, as the assembly 
process proceeds. 

The next step in the process is the incorporation of the insert 
equations in the system. As was described in detail in Section 3.3.3, the 
insert assembly consists of an insert by insert matrix multiplication and 
backsubstitution. Hie backsubstitution minimizes the number of equations 
required in the system since the displacement about the insert is 
eliminated from the system and only tractions are retained. 

A key problem in the entire process is the proper definition of 
appropriate coordinate systems, on a nodal basis. This is a problem common 
to any direct boundary element method which treats structures with 
nonsmooth surfaces. It arises because the tractions at a point are not 
uniquely determined unless the normal direction to the surface varies 
continuously at the point in question ♦ 

The original surface integral calculations are all done in global 
coordinates. If a displacement boundary condition is specified at a given 
node, in global coordinates, then no new coordinate system definition is 
required. It is only necessary to keep track of the subset of elements, 
containing the given node, on which the fixed displacement is to be 
reacted. However, if a displacement is specified in a nonglobal direction 
at a given node, then a new nodal coordinate system must be defined and, 
potentially, updated as further boundary conditions are processed. The 
associated nonzero reactions must then be expressed in the new coordinate 
system. 

Following this preparatory work, the final assembly of the system 
equations is carried out. It is performed in three major steps: 

1. Transformation of the columns of the matrices to appropriate local 
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coordinate systems and incorporation of any boundary conditions 
involving springs. 

2. Incorporation of compatibility and equilibrium conditions on 
interfaces between GMRs. On interfaces between two composite regions 
either a completely bonded condition (full displacement ccmpatibilty) 
or a sliding condition (only normal displacement compatibility) is 
available. At the interface between the insert and matrix a fully 
bonded connection has been assumed which will be relaxed later. 

3. Application of specified displacements and tractions. 

Two particular features of the equation assembly deserve special 
comment. First, in multi-GMR problems the system matrix is not full. 
Rather, it can be thought of as consisting of an NxN array of submatrices, 
each of which is either fully populated or completely zero. Only the 
nonzero portions of the system equations are preserved during system matrix 
assembly. In order to improve the numerical conditioning of the system 
matrix for the solution process, the columns are reordered to number 
variable lying on the same interface, but belonging to two different GMRs, 
as close together as possible. The rows of the system matrix are placed in 
the same order as the columns. 

Second, rather than simply assembling an explicit load vector at each 
time point in the solution process, load vector coefficient matrices are 
assembled and stored. These allow the updating of the load vector at any 
required time point simply by interpolating the time dependent boundary 
conditions and performing a matrix multiplication. A similar process is 
used in the calculation of interior and boundary stresses. 
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4.6 System Equation Solution 

The solver employed in Composite-BEST operates at the submatrix level, 
using software from the LINPACK package (Dongarra, 1979) to carry out all 
operations on submatrices. The system matrix is stored, fcy submatrices, on 
a direct access file. The decomposition process is a Gaussian reduction to 
upper trianglular (submatrix) form. The row operations required during the 
decomposition are stored in the space originally occupied by the lower 
triangle of the system matrix. Pivoting of rows within diagonal 
submatrices is permitted. 

The calculation of the solution vector is carried out by a separate 
subroutine, using the decomposed form of the system matrix from the direct 
access file. The process of repeated solution, required for problems with 
multi-time steps, is highly efficient. 
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5. PROGRAM DATA INPUT AID RESULTS 
5.1 Input Description 

The input to Composite-BEST is presently divided into five sections as 
follows : 

1. Case control (**CASE control) 

2. Material properites <**MATErial property) 

3. Generic modeling regions (**GMRegion) 

4. Interfaces (**INTErface) 

5. Boundary condition sets (**BCSEt) 

A detailed description of each of these sections is provided in the 
following paragraphs. The interface sections are optional; the other 

sections must be input at least once. 

Input quantities may be either alphanumeric or numeric (integer, 
floating point. E, or D format) as specified and may be up to 16 
characters. Individual entries on a card (both keyword and input) must be 
separated by at least one blank space. Input for certain keywords (as 
noted) may be continued onto more than one card by repeating the keyword on 

the new card(s) . 

Keywords may be input as shown; minimum input is the the CAPITALIZED 
characters. Those keywords which are underscored must always be input. 
Keywords shown below are indented to indicate groups of cards to be input 
together. However, it is not necessary to indent in this manner. 

The current program limits include: 

20 time points 

15 generic modeling regions 

600 elements (300 elements in problems having interior points) 
2500 nodes (560 nodes per region) 
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- 100 inserts per generic modeling region 
500 inserts per problem 

- 1200 source points (600 source points in problems having interior 
points ) 

302 source points per region in a local coordinate system 
99 interface element pairs (total) 

350 interface node pairs (total) 

60 boundary condition sets with springs 
maximum element number of 9999 
maximum node number of 9999 
maximum of 24 entries per input card 

1. Case Control Input 

Type of 

Keyword Input Input 

♦♦ CASE control 

Alphanumeric Case title 

Alphanumeric READ qi WRITE 

Numeric Output time value (s) 

The analysis is assumed to be static, constant temperature, elastic, 
and time independent unless the appropriate optional keyword is input. The 
optional keywords need be included only if a particular option is to be 
turned on. 

The case title should have a maximum of 72 characters. 

Input on the TIMES card may be continued on more than one card 

2. Material Property Input 

'Itie material property input section must be repeated for each separate 
material. 


TITL e 

RESTart 

TIMES 
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Keyword 

Type of 
Input 

Input 

♦♦MATErial property 



m 

Alphanumeric 

Material name 

TEMPerature 

Numeric 

Temperature value (s) 

EMODulus 

Numeric 

Young's modulus value (s) 

POISson 

Numeric 

Poisson's ratio value 

The Young's modulus must be input 

in the same order as the temperature 

values. 



Input on the TEMPerature and EMODulus cards may be continued on more 

than one card. 



NOTE: The elastic modulus of the 

inserts are defined in the GMR input. 

The Poisson 

Ratio for inserts is taken to be the same as the 

Poisson Patio 

of the Composite matrix. 

3. Generic Modeling Region Input 


The generic modeling region section must be repeated for each region. 

Keyword 

Type of 
Input 

Input 

*»GMR 



m 

Alphanumeric 

Region name 

MAT 

Alphanumeric 

Material name 

TREFerence 

Numeric 

Reference temperature value 

PQINts 

Numeric 

Node number, coordinate values 
(x,y,z) 

SURFace 

Alphanumeric 

Surface name, (reference 
surface name) 

TYPE 

Alphanumeric 

LINE qi CHAD 
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ELEMe nts 



Numeric 

Element number, node numbers 

EQBMal 

Alphanumeric 

Element number, + qi - 

INSErt 

Numeric 

Elastic modulus of insert 

POINts 

Numeric 

Node number, coordinate value 
(x,y,z) 

TYPE 

Alphanumeric 

LINEar or GUAD 

ELEMent 

Numeric 

Element number, radius of 
insert element, node numbers 

INTErior 

Alphanumeric 


POINts 

Numeric 

Point number, coordinate 


values (x,y,z) for interior 
sampling points 

Hie SURFace input may be either of two forms*. 

a TYPE card and an ELEMents card to define element connectivity 

The TYPE designation in SURFace input specifies the traction or 
displacement variation on the element. A surface may contain only one TYPE 
card. Therefore, if mixed variation is required in a region, two surfaces 
must be defined. 

Surface elements must have either 6 or 7 (triangles) or 8 or 9 
(quadrilaterals) nodes. Element numbering is consecutive around the 
boundary. 

Insert elements must have either 2 (linear) or 3 (quadratic) nodes. 
Nodes referenced in element connectivity must be explicitly defined under 
'insert' POINt card and the points should not intersect the outer boundary 
or other insert elements. 

The sign associated with the defining element on the NORMal card 
should be plus (+) if the element is numbered in a counterclockwise 
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direction as seen from the outside of the model or minus ( ) if it is 
numbered in a clockwise direction. Disjoint boundaries must have multiple 
element/sign pairs on the NORMal card. 

The points which are input in the INTERior input are treated as 
"interior" points. These points may be either nodal points, other 
surface points, or true interior points. 


4. Interface Input 

The interface input describes the connection of surfaces (or elements 
or points) of one composite natrix region to another. Interfaces between 
the composite matrix and inserts dfi D2t have to be defined. 

The interface section must be repeated for each interface. 

Type of 


Keyword 

Input 

Input 

■IMTErface 



GMR 

Alphanumeric 

Region rame of first region 

SURFace 

Alphanumeric 

Surface name in first region 

ELEMents 

Numeric 

Element number(s) in first 
region 

POINts 

Numeric 

Node number (s) in first region 

GMR 

Alphanumeric 

Region name of second region 

SURFace 

Alphanumeric 

Surface name in second region 

ELEments 

Numeric 

Element number (s) in second 
region 

POINts 

Numeric 

Element number (s) in second 


region 


SLIDing 


ijjjg interface is assumed to have complete displacement compatibility 
unless a SLIDing card is input, in which case only normal displacement 
compatibility is assumed. 
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The ELEMents card and/or the POINts card are included in SURFace input 

only to designate a subset of that surface. 

Input on the ELEMents and POINts cards may be continued on more than 

one card. 

5. Boundary Condition Set Input 

The boundary condition set section must be repeated for each new 


boundary condition. 

Type of 

Input 

Keyword 

iLEUt 

♦*BCSEt 



m 

Alphanumeric 

Boundary condition set name 

GMR 

Alphanumeric 

Region name 

VALUe 



RELAtion 



SURFace 

Alphanumeric 

Surface name 

ELEMents 

Numeric 

Element number (s) 

POINts 

Numeric 

Node number (s) 

TIMES 

Numeric 

Input time value (s) 

L0CA1 



GMR 

Alphanumeric 

Region name 

SURFace 

Alphanumeric 

Surface name 

ELEMents 

Numeric 

Element number (s) 

POINts 

Numeric 

Node number (s) 

DISPLacement 

Numeric 

Component value 

SPLIst 

Numeric 

Source point value (s) ci ALL 
or SAME 

T 

Numeric 

Time point identifier, dis- 
placement value (s) 

RIGId 

Numeric 

Component value 

SPRIng 

Numeric 

Component value, spring value 
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TRACtion 



SPLIst 

Numeric 

Source point value (s) at ALL 
or SAME 

T 

Numeric 

Time point identifier, trac- 
tion value (s) 


The ELE Ments card and/or the POINts card are included in SURFace input 
only to designate a subset of that surface. 

The TIMES card must be omitted in a boundary condition set which 
contains a RIGId card. If the value(s) on the TIMES card differ from those 
values in the case control input, the output is calculatd by linear 
interpolation. In the case of time independence (i.e., the TIMES card is 
emitted) the time point identifier on the T card must be l (one). 

The LOCA1 card designates input in the outward normal direction. The 
component value on the DISPLacement, RIGId, SPRIng, or TRACtion card must 
be l (one). Care must be taken not to mix global and local coordinate 
systems on a particular element. Care must also be taken not to input 
conflicting components on a particular node in a particular element. 

The VALUe card should be included with the Displacement card, the 
RIGId card, or the TRACtion card. The RELAtion card should be included 
with the SPRIng card. 

Either RIGId input, ar SPRIng input, as TRACtion input must be 
included in a boundary condition set. This input set may be included up to 
three times (once for each component) in a boundary condition set. 
However, different boundary condition types may not be mixed in a boundary 
condition set. 

The SPLIst card indicates the order in which the values are to be 
input on the T cards. The input may be in either of three forms: 

nodal values 
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ALL to indicate that a single constant value is to be input 
SAME to use the previous source point list within the current 
boundary condition set (this option may not be used for the first 
source point list in the current boundary condition set). 

Input on the ELEMents, POINts, SPLIst, and TIMES cards may be 
continued on more than one card. Input on the T card may be continued on 
more than one card, including the time point identifier on each card. 
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5.2 Sample Input 

Ihe following pages contain a sample data input for Compos ite-BEST for 
a cube with five inserts in tension. Note the simplicity of the input 
section for inserts. Each insert may contain a number of linear or 
quadratic element s. The keyword ELEMents is used before the definition of 
each distinct (non-connected) insert. 


**CASE 

TITLE CUBE WITH INSERT UNDER TENSION 
TIMES 1. 

RESTART WRITE 

* * H A TE 

ID M A T 1 
TEMP 70. 0 
EMOD 100. 

POIS 0. 3 

* * G MR 

ID GMR 1 
MAT MAT1 
TREF 70. 0 
POINTS 


1 

. 0000 

. 0000 

. 0000 

2 

. 5000 

. 0000 

. 0000 

3 

1 . 0000 

. 0000 

. 0000 

4 

1 . 0000 

. 5000 

. 0000 

5 

1 . 0000 

1 . 0000 

. 0000 

6 

. 5000 

1 . 0000 

. 0000 

7 

. 0000 

1 . 0000 

. 0000 

e 

. 0000 

. 5000 

. 0000 

1 001 

. 0000 

. 0000 

. 5000 

1 002 

. 5000 

. 0000 

. 5000 

1 003 

1 . 0000 

. 0000 

. 5000 

1 004 

1 . 0000 

. 5000 

. 5000 

1 005 

1 . 0000 

1 . 0000 

. 5000 

1 006 

. 5000 

1 . 0000 

. 5000 

1 007 

. 0000 

1 . 0000 

. 5000 

1 00B 

. 0000 

. 5000 

. 5000 

2001 

. 0000 

. 0000 

1 . 0000 

2002 

. 5000 

. 0000 

1 . 0000 

2003 

1 . 0000 

. 0000 

1 . 0000 

2004 

1 . 0000 

. 5000 

1 . 0000 

2005 

1 . 0000 

1 . 0000 

1 . 0000 

2006 

. 5000 

1 . 0000 

1 . 0000 

2007 

. 0000 

1 . 0000 

1 . 0000 

2008 

. 0000 

. 5000 

1 . 0000 


SURFACE SURF11 
TYPE QUAD 
ELEMENTS 


101 

1 

2 

3 

1 003 

2003 

2002 

2001 

1 001 

1 002 

1 02 

3 

4 

5 

1 005 

2005 

2004 

2003 

1 003 

1 004 

1 03 

5 

6 

7 

1 007 

2007 

2006 

2005 

1 005 

1 006 

1 04 

7 

8 

1 

1 001 

2001 

2008 

2007 

1 007 

1 008 

1 

1 

2 

3 

4 

5 

6 

7 

8 


201 

2001 

2002 

2003 

2004 

2005 

2006 

2007 

2008 



NORMAL 201 * 
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INSERT 1000. 


POI NTS 

3001 

0. 5 

0. 5 0. 

3002 

0. 5 

0.5 0.125 

3003 

0. 5 

0. 5 0. 25 

3004 

0. 5 

0.5 0.5 

3005 

0. 5 

0. 5 0. 75 

3006 

0. 5 

0. 5 0. 875 

3007 

0. 5 

0.5 1.0 

TYPE QUAD 


ELEMENTS 


301 

. 1 

3001 3002 3003 

302 

. 1 

3003 3004 3005 

303 

. 1 

3005 3006 3007 

INTERIOR 

POINTS 

9002 . 

5 . 1 

0. 5 

9004 . 

5 . 2 

0. 5 

9006 . 

5 . 3 

0. 5 


9007 .5.35 0.5 

901 0 .5 .375 0.5 

* * BC SET 

I D BC 1 
G MR G MR 1 
SURFACE SURF11 
ELEMENTS 104 
DISP 1 

SPLIST ALL 
T 1 0.0 

* * B C S E T 

I 0 BC2 
G MR GMR 1 
SURFACE SURF 1 1 
ELEMENTS 104 
POINTS 1001 1007 

D I SP 3 

SPLIST 1001 1007 

T 1 0. 0. 

* * BCS ET 

I D BC3 
GMR GMR 1 
SURFACE SURF 1 1 
ELEMENTS 104 
POINTS 8 2008 
D I SP 2 

SPLIST 8 2 00 8 
T 1 0. 0. 

* * BCSET 

I D BC4 
GMR G MR 1 
SURFACE SURF11 
ELEMENTS 102 
TRAC 1 
SPLIST ALL 

T 1 100.0 
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5.3 Output Description 

The output from Composite-BEST is relatively straightforward. It 
consists of ten sections, as follows: 

1. Complete echo of the input data set. 

2. Summary of case control and material property input. 

3. Complete definition for each GMR, including all surface insert nodes, 
surface and insert elements. 

4. Complete summary for each interface and boundary condition set, 
including the elements and nodes affected. 

5. Boundary solution (on an element basis), including displacements and 
tractions at each node of each element. 

6. The resultant load on each element and on the entire GMR is calculated 
and printed. 

7. Solution for the displacements and tractions at the Insert/Matrix 
composite interface (on an element basis). 

8. Displacement, stress and strain on a nodal basis, at all surface 
nodes, for each GMR. 

9. Displacements at interior nodes. 

10. Stresses at interior nodes. 

A sample output is shown on the following pages for the data input 
that was given in Section 5.2. 
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5.4 Sanple Output 


ORIGINAL pfiQP 

OF $»OOR QUALITY 


**** CASE CONTROL INPOT **** 


JOB TITLE CUBE WITH IHSERT ONDER TENSION 

TIHES TOR SOL0TION: 1.00000000 


BODNDART RESTART : 1 

BOUNDARY INTEGRATION EPSILON: 
INTERIOR INTEGRATION EPSILON: 


. 00100000 
. ooiooooo 


**** MATERIAL INP0T **** 

MATERIAL NAME: HAT1 

ELASTIC 

POISSONS RATIO: . 3000 

TEMP ALPNA E 

. 70000E*02 . 00000E*00 . 10000E*03 


**** GMR INPHT **** 


REGION 1 

NAME GHR1 MATERIAL HAT1 


PEEEPENCE TEMPERATORE 70. 00 

INITIAL TEMPERATORE OF GMR 00 


NOTES 57 ELEMENTS 6 

SOURCE P0INT3 45 CELLS 0 

HUMBER OF INSERTS 1 


COORDINATE LIST 


NODE 

I 

Y 

1 

1 

0000 

. 0000 

. 0000 

2 

. 5000 

0000 

. 0000 

3 

1 . 0000 

0000 

. 0000 

4 

1 0000 

5000 

. 0000 

5 

1 . 0000 

1 0000 

. 0000 

6 

. 5000 

1 0000 

0000 

7 

. 0000 

I 0000 

0000 

9 

. 0000 

. 5000 

. 0000 

1001 

. 0000 

. 0000 

. 5000 

1 002 

. 5000 

. 0000 

. 5000 

1 003 

1 0000 

0000 

. 5000 

1004 

1 0000 

5000 

. 5000 

1 005 

1 . 0000 

1 0000 

. 5000 

1006 

. 5000 

1 0000 

. 5000 

1007 

. ooon 

1 0000 

5000 

1 00R 

. 0000 

. 5000 

5000 

2001 

. 0000 

. 0000 

1 . 0000 

2002 

. 5000 

. 0000 

1 0000 

200 3 

1 . 0000 

0000 

1 . 0000 


S0PEACF3 1 

HOLE ELEHENTS 0 

INSERT ELEHCNT3 3 


40 



2004 1 

2005 1 

0000 

0000 1. 

5000 

0000 

1 . 0000 
1 . 0000 

ORIGINAL PAGE IS 

2006 

200? 

5000 1. 

0000 1. 

0000 

0000 

1 . 0000 

1 . 0000 

OF POOR QUALITY 

2008 

0000 

5000 

1 . 0000 


loot 

5000 

5000 

. 0000 


3002 

5000 

5000 

. 1250 


3003 

5000 

5000 

. 2500 


3004 

5000 

5000 

. 5000 


3005 

5000 

5000 

. 7500 


3006 

5000 

5000 

8750 


3007 

5000 

5000 

1 0000 


900 2 

5000 

1000 

. 5000 


900 4 

5000 

2000 

. 5000 


9006 

5000 

3000 

. 5000 


9007 

5000 

3500 

. 5000 


901 0 

5000 

3750 

. 5000 


COOROI KITE 

LIST OF NODES GENERATED BT CPBEST 


HOPE 

I 

r 

z 


73001 

5000 

6000 

. 0250 


7 300 2 

5000 

6000 

. 1 250 


73003 

5000 

6000 

. 2500 


73004 

5000 

6000 

. 5000 


73005 

5000 

6000 

. 7500 


73006 

5000 

6000 

. 8750 


73007 

5000 

6000 

. 9750 


83001 

4 134 

4600 

. 0250 


83002 

4 1)4 

4500 

. 1250 


8 3003 

4134 

4500 

2500 


n 3004 

.4134 

4500 

5000 


83005 

.4134 

4 000 

. 7500 


8 3006 

.4134 

4 500 

0750 


83007 

.41)4 

4 500 

. 9 ^ 50 


93001 

5066 

4500 

. 0250 


93002 

5966 

4 500 

, 1250 


93001 

5966 

4 500 

. 2500 


93004 

5966 

1500 

5000 


93005 

5 9*6 

4 500 

7 500 


93006 

5066 

4 500 

07 50 


9)007 

5966 

4 500 

9750 


SURFACE S0RF1 1 

QUADRATIC fARIATION 


ELEMENT 

NODES 




101 

1 

2 

3 1003 2003 

2002 2001 1001 1002 

102 

3 

4 

5 1005 2006 

2004 2003 1003 1004 

10) 

5 

6 

7 100- 7007 

2006 2005 1005 1006 

10 4 

7 

0 

1 1001 2001 

2008 2007 1007 1008 

1 

1 

0 

7 6 5 

4 3 2 

201 

2001 

2002 

2003 2004 2005 

2006 2007 2008 

INSERT ELEMENTS 

QUADRATIC fARIATION 


ELASTIC 

MOOHLOS OF 

INSERTS 

1 OOOOOE *0 4 


ELEMENT 

RADI US 

MOOES 


( I NSERT 

1 ) 




301 

1000 

3001 

3002 3003 


302 

1000 

3003 

1004 3005 


303 

. 1 0 n 0 

3005 

3006 3007 



SOURCE 

2 

POINT 

3 

LIST 

4 

5 

670 1001 

1002 1003 

1004 

1005 

1 006 

1007 

1009 

2001 

2002 2003 2004 

i 2006 

2007 

2000 

7 3001 

93001 93001 73002 03002 

97002 73003 

03003 

93003 

73004 

03004 

a 3004 

73005 

P30O5 o’nn*; 73006 

■ 93006 

71007 

93007 

93007 
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ORIGINAL HA3£ & 

OF POOR QUALITY 


mm BOUNDARY CONDITION IHFUT **** 
OMR GMR1 SURFACE SURF1 1 

ELEMENT LIST 

mm 


SOURCE POINT LIST 

7 g 1 1001 2001 

COMPONENT 1 DISPLACEMENT 

DATA VALUES: 

. OOOOOE *00 00000 E ♦ 0 0 

, OOOOOE *00 


2008 2007 

INPOT 

OOOOOE *00 


1007 1008 


OOOOOE * 00 


. 000008*00 


. 00000 E *00 


. oooooe*oo 


OMR OMR 1 


BOUNDARY CONDITION IHPOT **** 
SURFACE S0RF11 


ELFMFNT LIST 
104 

SOURCE POINT LIST 
1001 1 00 7 

COMPONENT 3 DISPLACEMENT INPOT 

DATA VALUES: 

OOOOOE *00 OOOOOE *00 


boundary condition IHPOT **** 
OMR OMR 1 SURFACE SURFM 

ELEMENT LIST 
104 

SOURCE rOIHT LIST 
8 2008 

COMPONENT 2 DISPLACEMENT INPUT 

data values- 

OOOOOE* 00 OOOOOE *00 


OMR GHR1 


boundary condttioh INFUT **** 
SURFACE SURF11 


ELEMENT LIST 
102 


SOURCE POINT LIST ^ ^ ^ ^ ^ 100< 


COMPONENT 1 TRACTION INPUT 


DATA VALUES: 

1 OOOOE *03 . 1 00008*03 

. 1 0000 E *0 3 


1 OOOOE *01 . 1 00008*03 . 10000E*03 


. 1 00 00 E * 0 3 . 1 0 0 0 0 F. ♦ 0 


matrix decomposition - 

CONDITION HUMBER 


diagonal block 

. 1 3387E*04 


1 


. OOOOOE *00 


j 1 0000 E * 0 3 
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ORIGINAL PAGE JS 
OF POOR QUALITY 


JOB TITLE- CUBE BITS INSERT ONOER TENSION 

BOONDAR Y SOLOTIOH AT TIME ■ 1.0000 FOR REGION » GMR1 


ELEMENT 

NODE NO 

I-DI5PL. 

Y-DISPL. 

X-DISPL 

X-TRAC. 

T TRAC. 

1 TRAC. 

101 

1 

. OOOOOE+OO 

. 1 5593E+00 

. 1 3754E+00 

. OOOOOE+OO 

. OOOOOE+OO 

. OOOOOE+OO 

101 

2 

. 48021 E + 00 

. 1 2997E+00 

. 1 3419E + 00 

. OOOOOE+OO 

. OOOOOE+OO 

. OOOOOE+OO 

101 

3 

. 98866E+00 

. 1 709 4 E *00 

. 1 4960E+00 

. OOOOOE+OO 

. OOOOOE+OO 

. OOOOOE+OO 

101 

1003 

. 97531 E + 00 

. 1 77 4 7 E *00 

. 32046E-06 

. OOOOOE+OO 

. OOOOOE+OO 

. OOOOOE+OO 

101 

2003 

. 98866E+00 

. 17094E*00 

1 4960E + 00 

. OOOOOE+OO 

. 00000 F +00 

. OOOOOE+OO 

101 

2002 

. 48021 E+00 

. 1 2997E*00 

1 34I9E + 00 

. OOOOOE+OO 

. OOOOOE+OO 

OOOOOE+OO 

101 

2001 

. OOOOOE+OO 

1 55?3E*00 

t 3754 E +00 

. OOOOOE+OO 

. OOOOOE+OO 

. OOOOOE+OO 

101 

1001 

. 00000E *00 

, 1 5477E*00 

. OOOOOE+OO 

. OOOOOE+OO 

. OOOOOE+OO 

. OOOOOE+OO 

101 

1002 

. 47753E+00 

. 1 37 3 4 E + 00 

. 57255E-06 

. OOOOOE+OO 

. OOOOOE+OO 

. OOOOOE+OO 


102 

3 

90866E+OO 

. 1 7094E+00 

. 1 4 960 E + 00 

. 1 OOOOE + 03 

. OOOOOE+OO 

OOOOOE+OO 

102 

4 

. 93793E+00 

1 4762E-02 

. 1 3883E+00 

. 1 OOOOE+03 

. OOOOOE+OO 

. OOOOOE+OO 

102 

5 

. 98891 E + 00 

1 74 26E + 00 

. 1 4 8 3<J E + 00 

. 1 OOOOE+03 

. OOOOOE+OO 

. OOOOOE+OO 

102 

1005 

. 97555 E +00 

1 7932E+00 

. 1 5470E-06 

. 1 OOOOE+03 

. OOOOOE+OO 

. OOOOOE+OO 

102 

2005 

98891 E+00 

- 1 7426E+00 

1 4 8 39E+00 

. 10000E+03 

. OOOOOE+OO 

OOOOOE+OO 

102 

2004 

. 93793E+00 

1 4761 E - 0 2 

1 3 88 3E +00 

. 10000E+03 

. OOOOOE+OO 

OOOOOE+OO 

102 

2003 

98865E+00 

. 1 7094E + 00 

1 4960E + 00 

. 1 0000 E + 03 

. OOOOOE + OO 

OOOOOE+OO 

102 

1003 

. 9753 1 E + 00 

. 1 77 4 7E + 00 

. 32046E-06 

. 1 OOOOE+03 

. OOOOOE + OO 

OOOOOE+OO 

102 

1 004 

. 92675E+00 

721 62E-03 

. 45097E-06 

. 1 OOOOE+03 

. OOOOOE+OO 

OOOOOE+OO 


103 

5 

. 98891 E+00 

1 74 26E+00 

. 1 4 8 39 E + 00 

. OOOOOE + OO 

. OOOOOE + OO 

OOOOOE+OO 

103 

6 

. 48072E+00 

1 331 3E + 00 

. 1 3248E + 00 

. OOOOOE + OO 

. OOOOOE+OO 

OOOOOE+OO 

103 

7 

. OOOOOE+OO 

1 5675E + 00 

. 1 3 6 96 E + 00 

. OOOOOE + OO 

. OOOOOE + OO 

. OOOOOE + OO 

103 

1 007 

. OOOOOE+OO 

1 5487E+00 

. OOOOOE+OO 

. OOOOOE+OO 

. OOOOOE + OO 

OOOOOE+OO 

103 

2007 

. OOOOOE+OO 

1 5675 E + 00 

1 3696E + 00 

. OOOOOE + OO 

. OOOOOE + OO 

. OOOOOE+OO 

103 

2006 

. 48072E+00 

-.1331 3 E +00 

1 3 2 4 7 £ ♦ 00 

. OOOOOE+OO 

. OOOOOE + OO 

OOOOOE+OO 

103 

2005 

. 9 98° 1 E+00 

1 7426E+00 

1 4 0 39E + OO 

. OOOOOE+OO 

. OOOOOE+OO 

OOOOOE+OO 

103 

1005 

97565E+00 

17932E + 00 

. 15470E-06 

. OOOOOE+OO 

. OOOOOE + OO 

. OOOOOE+OO 

103 

1006 

. 47808E+00 

1 378 7 E + 00 

4 7 4 4 1 E-06 

. OOOOOE + OO 

. OOOOOE+OO 

. OOOOOE+OO 


104 

7 

. OOOOOE+OO 

- 1 5676E+00 

. 1 3 6 96 E + 00 

89036E + O2 

. OOOOOE+OO 

OOOOOE+OO 

104 

8 

. OOOOOE+OO 

OOOOOE+OO 

1 401 4E+00 

- 10162E+03 

18829E+00 

OOOOOE+OO 

104 

1 

. OOOOOE+OO 

. 1 5693E + 00 

. 1 3 75 4 E+ 00 

- 9951 3E+02 

. OOOOOE + OO 

. OOOOOE+OO 

104 

1 001 

OOOOOE+OO 

. 1 54‘ , 7£*oO 

, OOOOOE+OO 

8901OE+O2 

. OOOOOE + OO 

- 1 051 4E O' 

104 

2001 

. OOOOOE+OO 

. 1 5593E+00 

1 3 75 4 E +00 

8951 3E+02 

. OOOOOE + OO 

. OOOOOE+f+O 

104 

2008 

. OOOOOE+OO 

OOOOOE+OO 

1 401 4E + 00 

- 10162E+03 

. 18820E+OO 

. OOOOOE+OO 

104 

2007 

. oooon e + oo 

1 5675E+00 

1 3 6 96 E ♦ 00 

89B36E+02 

. OOOOOE+OO 

OOOOOE+OO 

104 

1007 

OOOOOE+OO 

1 5 4 07 E+OO 

. OOOOOE+OO 

90021 E +02 

. OOOOOE + OO 

- 9 4 0 “ Q E 04 

104 

1000 

. OOOOOE+OO 

. 09522E-O4 

. 50108E-06 

-. 10605E+03 

. OOOOOE + OO 

. OOOOOE+OO 



JOB TITLE 

COBE 8 1 T H 

INSERT ONTER 

TENSION 





BOUNOARY 

SOLUTION AT TIME 

» 1.0000 

FOR REGION 

* GMR1 



ELEMENT 

NODE NO. 

I-DISPL. 

T-DISPL. 

1 - DISPL. 

X-TRAC. 

T TRAC 

Z TRAC 

1 

1 

. OOOOOE+OO 

. 15593E + 00 

. 1 3 7 5 4 E + 00 

. OOOOOE+OO 

. OOOOOE+OO 

OOOOOE+OO 

1 

8 

. OOOOOE+OO 

. OOOOOE+OO 

.1401 4E+00 

. OOOOOE+OO 

. OOOOOE+OO 

. OOOOOE+OO 

1 

7 

OOOOOE+OO 

1 5675E + 00 

1 3696E+00 

. OOOOOE+OO 

. OOOOOE+OO 

. OOOOOE+OO 

1 

6 

45072E+00 

1331 3E + 00 

. 1 3248E + 00 

. OOOOOE+OO 

OOOOOE+OO 

. OOOOOE+OO 

1 

5 

98091E+OO 

1 7 4 26 E+00 

1 4 8 39E+00 

. OOOOOE+OO 

. OOOOOE+OO 

. OOOOOE+OO 

1 

4 

93793E+00 

1 4* r 6?E-02 

. 1 300 3 E + 00 

. OOOOOE+OO 

. OOOOOE+OO 

. OOOOOE+OO 

1 

3 

. 99 866 E + 00 

. 1 7094 E + OO 

, 1 4 OKOE + OO 

OOOOOE+OO 

. OOOOOE+OO 

OOOOOE+OO 

1 

2 

. 48021 E+00 

1 2 99 7 E ♦ 00 

1 3419E + 00 

. OOOOOE+OO 

. OOOOOE+OO 

OOOOOE+OO 


201 

2001 

OOOOOE+OO 

. t 55O3E+00 

1 3 75 4 E + 00 

. OOOOOE+OO 

. OOOOOE+OO 

OOOOOE+OO 

201 

2002 

. 48021 E+00 

. 1 2997E+00 

1 341 9E + 00 

. OOOOOE + OO 

. OOOOOE + OO 

OOOOOE+OO 

201 

2003 

98966E+00 

. 1 7094E+00 

M960E+00 

. OOOOOE+OO 

. OOOOOE + OO 

OOOOOE+OO 

201 

2004 

93 7 0?E+00 

1 4 7 6 1 E * 0 2 

1 3I03E+OO 

. OOOOOE+OO 

. OOOOOE + OO 

. OOOOOE+OO 

201 

2005 

. 90891 E+00 

-. 1 7426E + 00 

-. 1 4939E + 00 

. OOOOOE + OO 

. OOOOOE + OO 

OOOOOE+OO 

201 

2006 

. 4B072E+00 

1 331 3E + 00 

1 3 2 4 7 E + 00 

. OOOOOE+OO 

. OOOOOE+OO 

. OOOOOE+OO 

201 

200' 

. OOOOOE+OO 

-. 1 5 67 5 E + 00 

- 1 3696E+00 

. OOOOOE + OO 

. OOOOOE+OO 

. OOOOOE+OO 

201 

2008 

. OOOOOE+OO 

. OOOOOE+OO 

1 481 4E + 00 

. OOOOOE+OO 

. OOOOOE + OO 

. OOOOOE + OO 


43 



JOB TITLE: CUBE HIT* INSERT UNDER TENS I ON 

load CALCULATION at TIME « 1.000000 


OKiCtffcAL ?A3E ?S 

OF POOR QUALITY 


LOADS FOR REGION GHR1 


ELEMENT 

X 

Y 

Z 

101 

. OOOOOE+OO 

. OOOOOE *00 

. OOOOOE*OQ 

102 

. 1 0000 E + 0 3 

. OOOOOE+OO 

. OOOOOE*00 

103 

. 00000 E *00 

. OOOOOE *00 

. OOOOOE+OO 

104 

99659E*02 

. 41 8I2E-01 

221 36E-04 

1 

. 00000 E *00 

. OOOOOE *00 

. OOOOOE+OO 

201 

. OOOOOE ♦ 00 

. OOOOOE+OO 

. OOOOOE+OO 

LOAD BALANCE 

. 31092E+00 

41 B42E-01 

- 22136E-04 

JOB TITLE: CUBE RTTR INSERT UNDER 

INSERT ELEMENT SOLUTION AT TIME « 

TEN3I0H 
1 . 0000 

FOR REGION * 


ELEMENT NODE HO. 


X-TRAC. T TRAC Z TRAC. 


301 

73001 

P9532E+Q0 

. 1 00 3 0 E < 01 

301 

83001 

. 1 2852E + 03 

-. 69704E+01 

301 

9 3001 

. 1 2607E + 03 

73593E+01 

301 

73002 

. ’2226E+00 

. 96271 E+01 

301 

8 3002 

. 1 1 060E+O3 

28546E+01 

301 

9 3002 

. 1 1989E*03 

29901 E»01 

301 

73003 

. 29744E + 00 

. 4 006 5 E +0 1 

301 

83003 

. 1 1 074E+O3 

171«7E*01 

301 

93003 

. 1 1 87 4 E *0 3 

2081 4 E+01 


302 

73003 

. 29 7 44E+00 

. 40O65E+O1 

302 

83003 

. 1 1 87 4 E+03 

171 67 E *0 1 

302 

9 3003 

. 1 1 074E + O3 

2001 4E+01 

302 

’3004 

. 1 4592E + 00 

. 64886E+01 

302 

83004 

. 1 1796E+03 

2 4 30 6 E * 0 1 

302 

93004 

. 1 1 74«E*03 

- 28693E+01 

302 

73005 

. 29 ’06 E + 00 

. 40O68E+O1 

302 

83005 

. 1 1 074E+OJ 

1 71 66E *01 

302 

93005 

. 1 1 874E*03 

2081 9E*01 


303 

73005 

. 29706E *00 

. 40O60E+O1 

303 

0 3005 

. 1 1 974E+0 3 

17166E+01 

303 

93005 

. 1 1 87 4E*03 

-. 2001 9E+01 

303 

73006 

. 72447E+00 

. 96234E+01 

30) 

03006 

. 1 1068E + O3 

2 05 3 4 E + 0 1 

303 

93006 

. 1 1 999E+03 

2 995 6 E ♦ 0 1 

303 

’3007 

. R°906E*00 

. 1 0029E *01 

303 

83007 

. 1 28 5 2 E *0 3 

- 69699E *01 

303 

9 3007 

. 1 260’E*03 

73594E + 01 


454I9E*02 
49540E+O2 
50428E+02 
10879E+02 
77024E+O1 
8 3 7 25E ' 0 1 
74677E+01 

- 88572E +01 

- B2263E+01 


- 746’7E«01 
-. B8572E+01 
-. 82263E+01 
. 27555E-03 
. 22600E '0 3 
. J6406E-03 
. 74674E+01 
. 805 75 E *0 1 
. 82249E+01 


. 74674E+01 
. 00 575E *01 
. 822 49E + 01 
. 1 0P85E*02 
. 77065E+01 
. 8301 6E+01 
. 4S43“’E*02 
. 4952‘ , E*0? 
. 5041 6E+02 


JOB TITLE: CM RE 0 I T H INSERT UNDER TENSION 

INTERIOR DISPLACEMENT AT TIME a 1.0000 FOR REGION * GMR1 


NODE 

I DISPLACEMENT 

Y DISPLACEMENT 

Z DISPLACEMENT 



«J002 

9001 

9006 

9007 
9010 

. 475566E+00 
. 4 7 30 ° 3 E + OO 
. 47079’E*00 
. 46«954E*00 
. 469622E+00 

. 1 05 7 7 7 E +00 
. ?2’4I’E-01 
. 362568E-01 
. 1 61 079E-01 
. 6761 20E * 0 2 

. 636808E-06 
. ’6599QE-06 
. 9 7<J ?25E~06 
. 1 1 4205E-05 
. 125234E-05 




JOB TITLE: CUBE 

INTERIOR STRESS AT 

RITR INSERT UNDER 
TIME * 1 0000 

TENSION 

FOR REGION 

* CMR1 


NODE 

SIGH A -IX 

SIGMA-YY 

SIGMA -ZZ 

TAO-XT 

TAU-XZ 

9002 

9004 

9006 

9007 
9010 

. 997037E+02 
. 1O1O03E+O3 
. 106541 E*03 
. 1 1 265BE+03 
. 116245E+03 

. 419311 E * 00 
. 1 31 69 9 E +00 
. 294053E-02 
04O79OE-O1 
35161 2 E + 00 

- 222506E+01 

3 4 34 3 1 E + 01 

4 220 3 5E*01 
451 955E+01 

-. 4 4 096 3 E + 0 1 

-. 074776E+O1 
032353E + O1 
725881 E+01 
636016E+01 
5501 59E+01 

627650E+01 
51 2052E*01 
375639E»01 
329232E+01 
323272E+01 


T AU-TZ 

6 7 540nE+O1 
624~66E*01 
52O628E+01 
424039E+O1 
-. 356346E+01 
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6. EXAMPLES OF COMPOSITE INSERT ANALYSIS 
6.1 Introduction 

In this section a number of examples are presented to verify and 
demonstrate the applications of the ceramic composite formulation for 
elastic fully-bonded inserts. 

In the mesh diagrams of the models containing the inserts, a double 
line is used to indicate the centerline of the insert elements. The length 
of these segments are shown in proper proportion for the three-dimensional 
views, however, the radii of the inserts are not indicated on these 
diagrams. The double line is a symbolic representation of the insert 
elements and does not in any way indicate the diameter of the insert. 
Refer to the example description for the values of the radii. 

Throughout this section consistent units are used in the definition of 
the examples. This means all lengths are defined in the same units and the 
tractions and the elastic moduli are defined in terms of these lengths as 
Force/ length 2 . No confusion should arise since the results are reported as 
non-dimensional quantities. 

6.2 Cube With a Single Insert 

The first test of the formulation is on a unit cube with a single 
insert through its center of radius 0.1. The cube is subjected to tension 
and shear in the direction parallel and perpendicular to the insert. The 
cube has a modulus of 100.0 and a Poisson ratio of 0.3. Consistent units 
are used for all information described in this problem. An insert with two 
different moduli of 1,000 and 10,000 is studied. The Poisson ratio of the 
insert is assumed to be the same as that of the cube. 

The problem is analyzed by both the present formulation and by a full 
three-dimensional multi- region BEM approach. As shown in Fig. 6.2.1, the 
model for the insert formulation consists of fourteen quadratic boundary 
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elements and the insert contains three quadratic insert elements. The two- 
region, three-dimensional model shown in Fig. 6.2.2 contain twenty 
quadratic boundary elements in the first region and sixteen in the second. 
Note 9-noded elements are used in describing the insert and hole to 
accurately capture the curvilinear geometry. 

In Fig. 6.2.3, the profile of the end displacement of the cube under a 
uniform normal traction of 100.0 (in parallel with the insert) is shown. 
The present formulation is in good agreement with the full three- 
dimensional results for Ei/E = 10. For the case Ei/E = 100, the insert 
formulation exhibits greater stiffness than the 3-D results. This 
difference is contributed by the way the load is distributed from the 
insert to the composite matrix. In the full 3-D model, the applied 
traction and the resulting reactions at the fixed end act directly on the 
end of the insert. In the composite formulation, the insert is assumed not 
to intersect the boundary surface and therefore the insert is moved back 
slightly from the end of the cube. The load is therefore transferred 
through the composite natrix to the end of the insert and to its sides in a 
manner that is slightly different from the full 3-D analysis. 

In Fig. 6.2.4, the stress distribution through the center of the cube 
(from A to B as indicated in the figure) is shown. Again the results are 
very good for Ei/E = 10, and deviates slightly from the full 3-D results in 
the second case. 

In Figs 6.2.5 and 6.2.6, the lateral displacements along the side of 
the cube are shown for a cube subjected to a shear traction of 100. For 
the case of applied shear perpendicular to the insert (Fig. 6.2.5), the 
results for both the insert and full 3-D model show good agreement. Once 
again a slight deviation is observed for Ei/E = 100. In the case of the 
shear traction in the plane of the insert (Fig. 6.2.6) the insert has 
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little effect on the displacement (as anticipated) and all results fall in 
close proximity. 

6.3 Lateral Behavior of a Cube With Multiple Inserts 

Existing methods of analysis of composite material based on mechanics 
of materials have been relatively successful in predicting the behavior of 
composite material for loading in the longitudinal direction. The 
properties perpendicular to the direction of the fibers are not so readily 
predictable by present means. The focus of the present example concerns 
this lateral behavior. 

Four cubes (Fig. 6.3.1) with one, two, five and nine inserts are fixed 
with a roller boundary condition on one side and subjected to a uniform 
traction, perpendicular to the inserts. The material properties, given in 
consistent units, are 

E insert . 10000 . jjnatrix , 100> 

^ insert = 0 3 v matrix = 0>3 

For the cube with one and two inserts, the boundary mesh consists of 
two quadratic surface elements on each lateral side and four elements on 
the top and bottom. For the cubes with five and nine inserts, one 
additional element was added to the side with the applied traction and to 
the side with the roller boundary condition. The top and bottom faces 
contain six elements to match the pattern of the sides. In all cases, each 
insert contained three one-dimensional quadratic elements. 

The profile for the end displacement for a cube with one insert and 
five inserts are shown in Figs. 6.3.2 and 6.3.3. The results seem to be in 
good agreement with the two-dimensional results. The 2-D results are 
approximations since plane stress is assumed. Hie 3-D solutions for the 
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one insert are within 2% error of the 2-D solution and within 3% for the 
case of 5 inserts. 

Also shown in Fig. 6.3.4 are the average end displacements for the 
one. two, five and nine inserts. Results show good agreement with 2-D 
results. For one, two and five inserts, the solutions are within 2% error 
of the 2-D results and 6% for the case of nine inserts where the insert of 
volume to total volume ratio is 28.2%. The result is also displayed in a 
plot of Effective Modulus vs. Insert Volume Ratio in Fig. 6.2.5. The 
effective modulus is defined as the average stress/average strain. The 
three-dimensional results follow closely to the two-dimensional solution. 

6.4 Thick Cylinder With Ciromferential Insert Supports 

The strength of a cylinder under internal pressure can be increased by 
adding stiff circumferential insert supports. In the present example, a 
three-dimensional, open ended thick cylinder with four inserts is analyzed. 
The inner and outer radii of the cylinder are 10 and 20 respectively, the 
thickness is 2 and the radius of the fully-bonded inserts is 0.5. By using 
roller boundary conditions on the faces of syirmetry, only a fifteen degree 
slice of the thick cylinder needs to be modeled. As shown in Fig. 6.4.1, 
sixteen eight-noded quadratic boundary elements are used to define the 
sides of the model, a nine-noded element is used on both the internal and 
external faces of the cylinder, and three insert elements are used per 
insert. Note, the inserts in this problem are curvilinear in geometry. 
The elastic modulus of the cylinder is assumed to be 100, and the effect of 
inserts with five different moduli of 100, 250, 500, 750 and 1000 is 
studied. The Poisson ratio is 0.3 for both the composite matrix and insert, 
and the internal pressure in the cylinder is 100. 

Results from a multi-region, axisymmetric BEM analysis were used for 
comparison with the 3-D insert results of the present example. The 


48 



axisymmetric model consists of twenty quadratic boundary elements on the 
outer surface, and six boundary elements per hole and per insert (Fig. 
6.4.2). The radial displacement through the thick cylinder along the top 
face is shown in Fig. 6.4.3 for all five moduli. Hie displacement for the 
composites with low Ei/E ratios are in good agreement with the axisymmetric 
results, and diverge slightly for higher Ei/E ratios. In Fig. 6.4.4, the 
circumferential stress is shown for the same points along the top edge. 
This stress is smooth for the homogeneous case (Ei/E = 1.0) and exhibits 
increasing fluctuations as the Ei/E ratio increases and the inserts take on 
more of the load. The circumferential stress of the 3-D insert model is in 
good agreement with the axisymmetric results for all cases. In Fig. 6.4.5, 
the radial stress is displayed for the two models. Hie inserts have little 
effect on this stress and the curves for the five moduli fall close 
together for both approaches. 

6.5 Cube With Multiple Inserts With Random Orientation 

In an attempt to analyze a material with a random fiber structure, 
cubes with multiple inserts oriented in random directions are studied. Hie 
cubes are of unit length and have four boundary elements per side (Fig. 
6.5.1a). Randomly oriented fibers of variable length with radii of 0.05 
are placed in five cubes in quantities of 5, 10, 15, 20, and 25 (Fig. 
6.5.1b-f). Three cases of material properties are considered for each 
cube. The modulus of the composite matrix is 100 for all cases, however, 
the modulus of the inserts are 500, 10,000 and 200,000 in the three cases 
studied. Poisson's ratio is uniformly 0.3 throughout. Roller boundary 
conditions are employed on three adjacent sides and a uniform normal 
traction of 100 is applied to a fourth face. 

The normal end displacement at the center of the face on the side with 
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the applied traction is plotted against the number of inserts in a cube for 
the three materials (Fig. 6.5.2). The displacement decreases with 
increasing number of inserts per cube and increasing Ei/E values as 
expected. 

6.6 A Beam With Insert Reinforcement in Bending 

In the last example, the applicability of the present formulation to 
the study of the micromechanical behavior of the ceramic composite is 
apparent. The present formulation, however, is equally applicable to 
typical problems encountered by civil engineers. Using Composite -BEST 
reinforced concrete can now be modeled exactly as a three-dimensional body 
and studied in detail for the first time. The present example considers a 
reinforced concrete beam. Here the concrete plays the role of the 
composite matrix and the reinforcement bars play the role of the fiber 
insert. In Fig. 6.6.1, a 4xlxl beam with four inserts is modeled using 
twenty-eight quadratic boundary elements. The ratio of insert modulus to 
matrix modulus (Ei/E) is studied for a range of valves between l and 100. 
The Poisson ratio is 0.3 for both the beam and reinforced rods. 

The beam is completely fixed at one end and a downward shear traction 
of 100 is applied to the other end. The non-dimensional vertical 
displacement of the end obtained from the analysis is shown in Fig. 6.6.2 
as a function of Ei/E. The non-dimensional displacement is defined here as 
the end displacement of the reinforced beam divided by the displacement of 
a homogeneous beam under similar conditions. 

The end displacement obtained from the mechanics of material solution 
is also displayed in Fig. 6.6.2 in non-dimensional form. The curvature of 
the two plots are very similar but differ in magnitude. This difference is 
contributed to the fact that although the mechanics of material solution 
accounts for the stiffening due to the inserts, it does not include the 
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effect of interaction between inserts. 


6.7 Laminated Fiber Composite 

A laminated composite fabricated from a fiber composite material is 
shown in Fig. 6.7.1. The fiber composite is constructed with a single row 
of fully-bonded fibers oriented in the same direction. A two-ply laminate 
is then constructed from the fiber composite with the fibers of the two 
layers oriented at 90° angles. A boundary element model created for the 
study of this material is shown in Fig. 6.7.2. A small slice containing 
two inserts in each layer is used. The model consists of two regions. The 
outer surface of each region is modeled with sixteen quadratic boundary 
elements and each insert contains two quadratic insert elements. The 
interface between the two regions is assumed to be a perfect bond, however, 
the present version of the program allows for sliding and spring 
connections also. 

The composite structure is subjected to bi-axial tension. This is 
accomplished with normal tractions of 100 applied to two adjacent roller 
boundary conditions applied to the opposite ends. The elastic modulus of 
the composite natrix of both regions are assumed to be 100, and the moduli 
of the inserts vary between 100 and 10,000. The Poisson ratio is 0.3 for 
both the composite matrix and inserts at all times. 

Figure 6.7.3 displays the displacement as a function of insert moduli 
for a point on the interface at the corner of the plate adjacent to the 
sides with the applied traction. The material exhibits less displacement 
as the modulus is increased, as expected. 
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Fig. 6.2.1 Discretization of an Insert in a Unit Cube Utilizing 
Quadratic Insert Elements 



Fig. 6.2.2 Full Three-dimensional, Multi-region Discretization 
of an Insert in a Unit Cube 
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Fig. 6.2.3 Comparison of Displacement Profiles Between the Full 
3-D Model and the Insert Element Model for the End of 
a Cube in Tension 
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Fig. 6.2.4 Axial Stress Through the Cross Section of a Unit Cube 
in Tension with a Single Insert in Parallel with the 
Loading 
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Fig. 6.4.3 Radial Displacement Through a Pressurized Thick 

Cylinder with Circumferential Inserts for Ei/E = 1.0, 
2.5, 5.0, 7.5, 10.0 



Fig. 6.4.4 Circumferential Stress Through a Pressurized Thick 

Cylinder with Circumferential Inserts for Ei/E = 1.0, 
2.5, 5.0, 7.5, 10.0 
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RADIAL STRESS / INTERNAL PRESSURE 



1.0 1.5 2.0 


RADIAL DISTANCE / II^NER RADIUS 


Fig. 6.4.5 Radial Stress Through a Pressurized Thick Cylinder 
with Circumferential Inserts for Ei/E = 1.0, 2 5, 
5.0, 7.5, 10.0 
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6.5.2 End Displacement of a Unit Cube With Random Oriented 
Inserts of Ei/E = 5.0, 10.0 and 100.0 


Fig. 6.6.1 Discretization of a Reinforced Beam Utilizing 

Quadratic Insert Elements to Model the Four Inserts 
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Fig. 6.6.2 Non-dimensionalized Vertical End Displacement of a 
Reinforced Beam in Bending vs. the Modulus of the 
Insert Over Modulus of the Beam 



Fig. 6.7.1 Laminate-Fiber Composite 
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7. SUMMARY OF CURRENT ACHIEVEMENT 

As is evident from the previous section, significant progress has been 
made towards the goal of developing a general purpose boundary element 
program for the micromechanical studies of advanced ceramic composites. 
Ihe formulation for elastostatic analysis of fully-bonded inserts has been 
implemented and validation runs have been shown to correlate with results 
obtained by full three-dimensional boundary element analyses. The 
composite insert code development is based on the advanced boundary element 
program BEST3D and all of its general purpose features have been retained 
in the new code. These facilities which include definition of local 
boundary conditions and multi-region substructuring will allow real 
problems encountered in industry to be analyzed. 

Considerable effort has been focused in developing a formulation that 
is not only accurate, but computer efficient. For instance, analytical 
integration has been performed on the kernel function of the insert and 
hole in order to expedite numerical integration. Also a new boundary 
integral equation formulation was developed to both facilitate an efficient 
assembly scheme for the inserts, and to reduce the number of unknowns in 
the system and therefore render a smaller set of equations which is less 
expensive to solve. 

Overall the cost of the new composite insert analysis is just slightly 
more expensive (for a moderate number of inserts) than the cost of 
analyzing the matrix without the inserts present. The additional cost is 
primarily attributed to the additional integration of the outer surface of 
the matrix which is required for additional nodes on the insert's hole, and 
towards the expense of solving a slightly larger system. In any case the 
price is far less than the cost required for a full three-dimensional, 
multi-region analysis of the same problem since in the 3-D approach more 
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two-dimensional surfaces must be integrated and far more additional 
equations (for the nodes describing the hole ao£ the insert region) must be 
integrated and added to the equation system. The data preparation for the 
present method is far less involved than a full three-dimensional modeling 
of the composite insert using the ordinary multi-region approach. 
Furthermore, the code developed for the present work allows up to 500 (100 
per GMR) insert elements in an analysis. An ordinary multi-region code 
would require a 500 region capability in order to compete. In terms of 
computer expense and the cost of da t a preparation for a 500 region problem, 
such analyses would be impractical. 

The work in this initial year has demonstrated the accuracy and 
efficiency of the composite insert formulation applied to elastostatics. 
The method is not only perhaps the best tool, but also may be the only 
practical tool for the analysis of composite inserts in real problems 
encountered in industry. The boundary element method already has been 
proven successful in multi-region, transient and steady-state, 
elastodynamic and heat conduction analysis. Goupled with the success of 
the present work, the plan of this contract to extend the composite insert 
formulation in these other areas holds great potential. 
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8. FUTURE DEVELOPMENT 

Significant progress has been achieved in the first year in 
developing a user friendly# ceramic composite base program. The 
infrastructure of the program is set up to facilitate future development. 
Presently the elastoplastic analysis for fully-bonded composite inserts has 
been implemented and tested. Delivery of the present code is planned for 
March 1989. 

The remaining work on the present contract consists of: extending the 

formulation to steady-state and transient elastodynamic and heat transfer 
analysis; developing sophisticated interface connections including failure 
at the interface; and incorporating thermal and nonlinear material 
phenomena in the problem. The sequence for development of these tasks are 
inconsequential since they are relatively independent of one another. A 
workplan for their development is presented belcw, however, coirmunication 
with interested parties at NASA will influence the exact order in which the 
work will be performed. 

The primary thrust for the upcoming year will be directed toward the 
incorporation of more sophisticated interface conditions. Phenomena of 
interest include imperfect bonding, progressive debonding, cracking of the 
matrix, and controlled tension failure of the fibers. Initially this work 
will be developed for elastostatics, and extended to the other analysis 
types as they are developed. 

Also during the next year, work will focus on the development of a 
ceramic composite analysis capability for steady-state heat conduction and 
the extension of thermal effects in elastostatics analysis. One- 
dimensional representations of fully-bonded ceramic fibers will be 
developed and incorporated into the code for heat conduction analysis, 
however, the suitability of the one-dimensional functional approximations 


67 



roust be critically examined for geometries and material properties typical 
of ceramic composites. In particular, it may prove necessary to employ a 
sine-cosine variation in the circumferential direction similar to the 
elastostatic formulation. Accuracy will be assessed via comparison with 
full three-dimensional BEST3D analyses. During the development for steady- 
state heat conduction and thermoelasticity, the thermal contact resistance 
between the matrix and insert interface will be added. 

Presently the insert fibers are entered in the data input of the 
program as individual one-dimension line element with a prescribed radius. 
This provided tremendous advantages over the full 3D approach in terms of 
both modeling effort and computational efficiency. However, from a user's 
standpoint, even the one-dimensional representation is cumbersome if each 
fiber centerline must be defined. Clearly, a more convenient user 
interface is needed in order to provide a practical tool for 
micromechanical analyses of ceramic composites. This interface has been 
designed and will be implemented as part of the second year effort. As a 
result, the positioning of fibers will be determined internally by the 
program within any generic modeling region based upon the specification of 
a few keyword-driven parameters. With this approach, the user need only 
define the outer surface of the body, as in any typical BEST3D model, along 
with parameters specifying the overall fiber content, shape, size and 
orientation. Finally, an updated version of the ceramic composite code 
will then be made available at year end. 

During the remaining three years, ceramic composite analysis 
capability will be developed for transient heat conduction, steady-state 
and transient elastodynamics, and material nonlinearities at high 
temperatures. Ihe primary task will involve the development of appropriate 
representations for the inserts. Of course, now this task is much more 
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difficult due to the complexity of the governing differential equations and 
fundamental solutions. Additionally, investigations will be required to 
determine appropriate nonlinear material models and solution algorithms for 
the elevated temperature response of ceramic composites, and extensive 
testing of all of the capabilities will be conducted. Specific priorities 
for this advanced work will be established during 1989. Validation and 
verification runs will be conducted regularly, so that reliability is 
maintained throughout the duration of the five-year program. Once again, 
at the end of each program year, the resulting general purpose 
micromechanical ceramic composite code will be deliverable. 

The final version, in particular, will provide a very precise, yet 
very efficient, user-friendly, design and analysis tool for ceramic 
composites exposed to severe operating environments. The resulting 
computer code will enable an engineer to undertake rapid numerical 
experiments to gain insight into the micromechanics of a particular 
composite. Armed with this information, the disposition of fibers can be 
selected to optimize performance under inelastic, thermal and dynamic 
loading. 
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LIST OP SYM3CLS 




C ij(U 


G ij' F ij 


N y (h) 




u i 


X 


X 


y 





v 


5 


ij 


Boundary system matrices in assembled form 

A tensor dependent on location of the field point ? 

Kernels of the displacement equation 

One-dimensional shape function 

Two-dimensional shape function 

Traction 

Displacement 

Refers to global coordinates of an integration point 
System vector of unknown boundary quantities 
System vector of known boundary quantities 
Strain 

Refers to local coordinate of an integration point 
Poisson's ratio 

Refers to coordinates of a field point 
Stress 


Subscript 

, Spatial derivative 


i.j.k 


Indie ial notation 
i.j.k = 1.2,3 
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